Thermalization of Strongly Disordered Nonlinear Chains 
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Thermalization of systems described by tiie discrete non-linear Schrodinger equation, in the strong 
disorder limit, is investigated both theoretically and numerically. We show that introducing corre- 
lations in the disorder potential, while keeping the "effective" disorder fixed (as measured by the 
localization properties of wavepacket dynamics), strongly facilitate the thermalization process and 
lead to a standard grand canonical distribution of the probability norms associated to each site. 
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Introduction - The discrete nonlinear Schrodinger 
equation (DNLSE) continues to be an active area of re- 
search. It became a paradigm in the field of nonlinear dy- 
namical systems pL^ and it occupies an important place in 
nonlinear optics [5] and the theory of ultra-cold atomic 
Bose- Einstein condensates [51. Much of recent work is 
dealing with the analysis of a wavepacket spreading, i.e. 
only one site (or, perhaps, a compact group of a few) is 
excited and the time evolution of this initial excitation 
is followed. The focus of this activity is the interplay 
of non-linearity and Anderson localization induced by a 
random potential jl| . A different type of problem jB] E] is 
when, initially, all sites are excited, so that the amount of 
energy endowed to the system is proportional to its (in- 
finite, in the thermodynamic limit) volume [7]. Follow- 
ing the evolution of such initial preparation one can ask 
whether a large isolated system thermalizes. Specifically, 
one can consider a large but finite part of the isolated 
system and ask whether it approaches a grand canonical 
distribution with some values of the temperature T and 
chemical potential fi. 

The two types of problems, i.e. "spreading" vs "ther- 
malization", are physically distinct, as emphasized in the 
work of Basko j5]. In the former case, the system is ex- 
cited locally (a microscopic excitation), whereas in the 
latter the excitation is global i.e. macroscopic. In a fi- 
nite size system the distinction between the two prob- 
lems can become blurred, if one does not specify how the 
system's energy scales with its size. It is possible, for 
instance, to consider the case when the energy provided 
to the system remains fixed and finite, while the size of 
the system increases [SI [S] (the problem of a wave packet 
spreading belongs to this category). In such a case the 
system, assuming it thermalizes, will reach a tempera- 
ture which is inversely proportional to system's size and, 
thus, approaches zero in the thermodynamic limit. 

Here we consider genuine, "thermodynamic" thermal- 
ization, when the system's energy is proportional to its 
size. We concentrate on the case of strong disorder, 
when in the absence of the non-linearity, all modes are 
strongly localized. We consider both uncorrelated and 
correlated random potentials and show that introduc- 



ing correlations in the disorder, while keeping the effec- 
tive disorder fixed (as it is measured by the spreading of 
an initially localized wavepacket in the absence of non- 
linearity), strongly facilitate the thermalization process. 
Model- The time-dependent DNLSE is given by 



dt 



^Vn - viljjn+l + -ipn-l) + XlV'nPV'n 
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where e„ is the n-th site energy, while the coupling con- 
stant V is set below to unity. The nonlinearity strength 
X and the time t are measured, respectively, in units of 
V and v~^ . The model ([l]) can be interpreted as a sys- 
tem of coupled nonlinear classical oscillators [3 [10] , with 
each oscillator being identified with the corresponding 
lattice site. We use "site" and "oscillator" interchange- 
ably, and the same goes for "site energy" and "oscilla- 
tor frequency". Below we will consider a finite lattice 
1 < n < N with periodic boundary conditions. To define 
the problem completely one has to specify the sequence 
of e„'s as well as the initial preparation. 

We will study three kinds of e-sequences: (A) The 
step-potential where the first N/2 sites are assigned on 
-site energy eL = 0, while the rest N/2 sites are hav- 
ing energy cr = const ^ 0; (B) Uncorrelated random 
potential where the e„'s are independent random vari- 
ables chosen from a rectangular distribution of width 
2W, centered at zero; and (C) Correlated random po- 
tential where a random sequence of en's, such as in (B), 
is created first and then it is convoluted according to: 
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where D 



X;fe[(l + cos(27rfc/7V))/2]^*^/'^'. This choice of D en- 
sures that the second moment of the correlated potential, 
(e^), is equal to that of the corresponding uncorrelated 
one, (e^). A simple calculation shows that, for N/l ^ 1, 
the correlation function (e„e„+m) ~ (e^) exp(— m^/C^), 
where the correlation length is C = ly/2/Tr. 

The initial wavefunction amplitudes {tpn} associated 
to the n— th site, were chosen for all cases to be a random 
number in the interval [0.75,1.25], while the normaliza- 
tion JV = J2n IV'nP = ^ ■^^s imposed. Note that this 
implies that, initially, the phases of all V'n's are zero. In 



addition to the conserved norm M, the total energy of the 
system H is also conserved by the dynamics of Eq. (fTl) : 

n = Y.[enHn\' - t'lV^^n+l + ^c) + Jxl^j'], (2) 



where the three terms on the r.h.s correspond, respec- 
tively, to potential, kinetic and interaction energy. The 
conserved quantities were monitored in our simulations 
to insure accuracy of the 4rth order Runge-Kutta scheme. 
Theory - We will investigate the thermalization pro- 
cess for the system of oscillators described by Eq. ^ . We 
assume that the coupling between oscillators is weak. For 
the random cases, (B,C), this is ensured by strong disor- 
der, i.e. by the large parameter W/v. For the step po- 
tential there is no randomness and the effective coupling 
strength depends on the initial preparation. It turns 
out that the aforementioned preparation corresponds to 
a high temperature regime, in which case the effective 
coupling between oscillators is weak [SI [?]• Under the 
weak coupling condition the thermalization criterion can 
be applied to an individual oscillator: in equilibrium, an 
oscillator with eigenfrequcncy e„ should obey the grand 
canonical distribution (GCD) for its norm /„ = l-fAnP, 
with the phase of ipn being completely random [5] : 




FIG. 1: (Color online) Distribution P{I) of the site-norm / 
for the case of a lattice with A^ = 100 sites, x — li ^nd a step- 
potential with eL = and er, = 7.5. Black (o)/red (•) colors 
correspond to left/right sites. The solid lines are best fits to 
Eq. (^. The fitted values {/Sl = 0.258; /it = -0.109) for 
the left and (/3_r = 0.253; /ijj = —0.1) for the right sub-chains 
respectively, agree nicely (within numerical accuracy of the 
fit) with one-another. Inset: Evolution of the first moment 
of the time-averaged on-site norms (/„(f)). The dashed lines 
are the theoretical predictions of Eq. (|6k) with /3 = 0.258 
and fi — —0.109. Blue lines indicate border sites n = 1 and 
n = 50. 
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Above /? is the inverse temperature, while the partition 
function Z„ is given by the integral 

Zn = / din exp[/3(^ - en)In - -/?x/^] 

^eM^WfcU^En). (4) 

2^X 2x y 2x 

In Eq. Q, we have introduced the notation £'„ = en — fJ,- 
Instead of reconstructing the full distribution, we will 
focus our investigation (except from the case (A)) at the 
time average of the norm and its second moment of each 
individual oscillator: 
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Their time dependence will allow us to conclude whether 
or not a specific oscillator has reached equilibrium. In 
the long time limit, these quantities should be equal to 
their statistical average, with respect to the GCD given 
by Eqs. (3J4|. The latter averaging results in: 
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For sufficiently weak nonlinearity, when x ^ PEn, the 
partition function Z„ can be expanded, leading to: 
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Numerical Analysis - We now discuss our numerical 
results for each of the three on-site potential sequences. 
We start with the case (A) associated with the step po- 
tential. This is an instructive example of thermalization. 
We recall that the chain of N sites is divided into two 
equal parts- the left and the right sub-chains with site 
energies e^, = and e^ ^ 0, respectively. The initial 
preparation with random on-site norms /„ ^ 1 evolve 
with time according to Eq (II]). Thermalization of the 
system is achieved due to a significant redistribution of 
the norm between various sites: "flow of the norm" from 
the right sub-chain (large site energies) towards the left 
sub-chain (zero site energies) occurs during the thermal- 
ization process. At equilibrium, a strong correlation be- 
tween the norm at a specific site n and the potential 
energy e„ is developed, so that the requirement Eq. ^ 
(or the more detailed condition for the full distribution) 
is satisfied. This scenario is confirmed by our numerical 
data for (/„)(i), reported in the inset of Fig. [Tlfor a chain 
of A^ = 100 sites with x = 1 and eright — 7.5. Indeed, for 
short times (roughly up to t = 100) the oscillators are 
exchanging norm among them and the various curves 
fluctuate wildly and cross each other. At larger times 
they " get organized" into two groups of curves (with the 
exception of the few border sites between the two sub- 



chains). Note, however, that for a long interval of time 
(from t ^ 100 till i '^ 5 x 10*) the two groups are very 
close to one another and do not change with time. For 
still longer times, they start diverging and eventually sat- 
urate at two different values for (/„) which are predicted 
by Eq. ^. We interpret this behavior as two-stage ther- 
malization. First, oscillators within each sub-chain ex- 
change norm and energy among them, with very little 
interaction between the two groups of oscillators. This 
rapid process (on time scales t ~ 100), leads to "internal" 
thermalization, associated to each sub-chain separately. 
Then, at much longer times, t ~ 5 x 10*, the slow pro- 
cess of interaction between the two groups of oscillators 
becomes efficient and it eventually leads to a full equilib- 
rium. Thus, the potential discontinuity, between the two 
sub-chains, serves as a "bottleneck" which slows down 
the total thermalization and provides the long-time scale 
for the thermalization process. 

We have also investigated the whole probability dis- 
tribution P{In) for the step potential chain. After an 
initial transient t ^ 10^, we have collected the norm /„ 
of each oscillator at various times. In all cases more than 
10^ data were used for statistical processing. We have 
found that the oscillators follow two distinct distribu- 
tions, which are characterized by their on-site energy e 
(i.e. left or right sub- chain) according to the prediction 
of Eq. (Is]). We have observed that the oscillators that 
correspond to high on-site potential energy were first to 
thermalize while the ones on the left sub-chain, corre- 
sponding to e„ = 0, achieve thermalization with a slower 
rate. A representative example of P{In) is shown in the 
main panel of Fig. [T] For better statistical processing we 
have averaged the left and right distributions of all oscil- 
lators. We have found that both distributions P^^ (/) and 
Pejj(/) follow nicely Eq. (pi. The corresponding fitted 
left and right chemical potentials and temperatures are 
in excellent agreement with one another, indicating that 
the total chain reached global thermalization. 

(B) Next, we investigate the thermalization of a chain 
with uncorrelated on-site potential. We concentrate on 
the analysis of the first (/„) and second moment (I,'^) 
of the norm. We consider cases for which the coupling 
parameter, r = W~^ and the effective non-linearity 
strength p = xIn/W are small 0. We let the initial 
preparation to evolve for a long time, of the order of 
5 X 10^ and then register the value of the norm (/„), for 
all N oscillators. Representative results of (/„) versus 
the site energies e„'s ior W = 6, N = 100, x = 1 a-^d two 
different initial preparations are shown in Fig. [2k,c for 
times t = 5x 10^ 7.5 x 10^ 10^; 2x 10^ 5x 10^ Such plot 
provides a check for the validity of Eq. (|6| which should 
hold in equilibrium. In Fig. [2]3,d we report also the sec- 
ond moment, {1^)- One can see that for high-energy sites 
^1 < £« < 6 both quantities, already for t = 5 x 10^, 
follow a regular curve and stay there for longer times. 
This part of the data is nicely fitted to the expression 
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FIG. 2: (Color online) First {I„{t)) (a,c) and second {In{t)} 
(b,d) moment of the on-site norm vs. e„ for an uncorrelated e- 
sequence with W — 6. We considered a lattice with A*' — 100 
and X = 1- Subfigures (a,b) and (c,d) correspond to two 
different initial preparations. Solid lines correspond to Eq. 
(|6k,b) 



Eq. ([6]), with appropriate /? and ^ values. Moreover, we 
have checked that for e > 4 the norm (/„) and, thus, the 
nonlinearity parameter become sufficiently small for the 
simplified relations Eq. (ItI) to hold. We have found, for 
instance, that for these highest energy tails, the relation 
(/^) = 2(/„)^ is satisfied quite accurately. This is an in- 
dication for a negative-exponential distribution (Eq.([3]) 
without the nonlinear term). 

We now turn to the low energy part, e < —3, of the 
data in Fig. [2] These sites do not show a clear sign 
of thermalization. The points do not fall on a regular 
curve but rather form a "blob" . Moreover, in the course 
of time the points are "jumping around" with no clear 
indication of approaching a regular curve. We conclude 
that, while the high energy sites thermalize already for 
relatively short times, the low energy ones do not ther- 
malize for much longer times (of course, they well might 
thermalize for even longer times). Let us note that for the 
lowest energy sites the nonlinearity parameter p reaches 
a value of about 1/2, so that the weak nonlinearity con- 
dition 5J is not really satisfied. 

(C) Finally, we discuss the case of a correlated disor- 
der and demonstrate that the thermalization process for 
this case is quite different from that in (B). A random 
correlated sequence of e„ is created as explained above 
(from now on we remove the "tilde"). The random cor- 
related ensemble is characterized by the "bare" disorder 
W and the correlation length (. The first thing to note 
is that it would not be appropriate to compare uncorre- 
lated and correlated ensembles with the same W; the 
point is that introducing correlation while keeping W 
fixed will, generally, change the "effective disorder" as 
measured, for instance, by spreading of an initially local- 



ized wave packet. Therefore the statement correlations in 
the randomness facilitate thermalization becomes mean- 
ingful only when the same effective disorder is kept for 
the cases (B) and (C). As a measure of effective disorder 
we chose the spreading of a particle (for x = 0) , initially 
placed at some site. More precisely, we consider the en- 
semble averaged (log jV'raP), in the long time limit, where 
n is the distance from the initial site. This quantity is 
plotted in Fig. |3]3 for two different ensembles: uncor- 
related disorder (case (B)) with W = 6 and correlated 
disorder (case (C) with W = 7 and ( = 2.7). Comparing 
the two curves one can conclude that, indeed, these two 
ensembles have the same effective disorder. 

We are now in the position to study thermalization in a 
chain with correlated disorder and to meaningfully com- 
pare the results with those in Fig. [2] for the uncorrelated 
case. In Fig. l3k,b we report our numerical results for the 
correlated disorder, in the same fashion and for the same 
successive times as in Fig. [2j Unlike Fig. [2j now all the 
data, including the low energy sites, fall on one smooth 
curve, given by Eq. (|6| with /? = 0.11 and /i = —10. Note 
that the spread in the site energies here is significantly 
larger than in Fig. [2J This happens because in the con- 
volution process e„'s larger (in absolute value) than the 
maximal "bare" value {W = 7, in our case) emerge. This 
wide spread of e^'s explains why the chemical potential 
is so low. Indeed, if the nonlinearity were very small, 
then Eqs. ([7| would imply that ^ must be considerably 
lower than the lowest site energy. Although in our case 
the nonlinearity is not very small, nevertheless, the large 
spread in the site energies pushes the chemical potential 
down (provided, that the system does thermalize). 

In Fig. |3j:,d we plot the numerical results for the same 
disorder realization as in Fig. [3k, b but with a different 
initial preparation. These data display an interesting fea- 
ture: the green rombs or the blue triangles arrange them- 
selves nicely in two distinct curves, which merge into one 
for larger e's. This suggests that we have here two sub- 
systems of sites which thermalized separately, with two 
somewhat different /3 and fj,. And then, in the course of 
time, there is a slow process of equilibration between the 
two subsystems. It is tempting to suggest that the points 
which form the upper or lower curves originate from sites 
belonging to one or the other half-chain. We have checked 
and confirmed this hypothesis. The asymmetry between 
the two half-chains is due the fluctuations in e's, as well 
as in the initial values of /„'s. We have also studied some 
other disorder realizations from the aforementioned cor- 
related ensemble. While the details of the time evolution 
depend on the realization and the initial preparation, we 
have observed in all examples an approach to equilibrium 
- in sharp contrast with the uncorrelated case. 

Conclusions - We demonstrated that correlations in 
disorder strongly affects thermalization. We can offer the 
following intuitive explanation for this effect: In the ab- 
sence of correlations (and for strong disorder) the mecha- 




FIG. 3: (Color online) Same as in Fig. [2] but for a corre- 
lated e„-sequence with W — 7 and correlation length ( = 2.7. 
Subfigures (a,b) and (c,d) correspond to two different initial 
preparations (the same as used for Fig. [2a,b and c,d respec- 
tively). In subfigure (e) we report the asymptotic wavefunc- 
tion for a linear lattice with an uncorrelated (black dashed 
line) and a correlated (red solid line) e -sequence with W = 6 
and W = 7 (and ( = 2.7) respectively (the initial condition 
was a single site excitation). 



nism of thermalization occurs via "resonant triples" , i.e. 
clusters of three oscillators with close eigenfrequencies e„ 
[5]. Correlated disorder enhances the probability of oc- 
currence of such resonant triplets, thus, facilitating ther- 
malization. We should stress again, however, that, since 
in our study the strong inequalities r << l,p << 1 are 
not satisfied, we are not quite in the regime of a very 
weak coupling and nonlinearity, studied in Ref. [5]. 
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